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The influence of local dephasing in the transition of ballistic to diffusive transport is investigated. 
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I. INTRODUCTION 

Ever since the discovery of anomalous heat transport 
though a chain of coupled harmonic oscillators in the 
seminal work of Rieder, Lebowitz, and Lieb [T], the topic 
of heat transport through systems of harmonic oscillators 
and the question of how classical diffusive heat transport 
emerges from a microscopic classical or quantum descrip- 
tion remains an interesting field of research. The classical 
law of diffusive heat conduction as stated by Fourier [2] , 

J = -kVT, (1) 

relates the heat current to the negative gradient of the 
temperature by the thermal conductivity k, the latter 
being a property of the material. Fourier's law of heat 
conduction contains several key statements: i) at con- 
stant system size, the magnitude of the heat current is 
proportional to the temperature difference (J oc AT), 
ii) at constant temperature difference and constant ther- 
mal conductivity, the heat current scales inversely pro- 
portional to the distance between the heat baths, which 
is given by the system dimension separating the baths 
(J oc l/L). Rieder, Lebowitz, and Lieb found for their 
one-dimensional classical system that the heat current is 
proportional to AT but does not scale with the system 
size. Nonetheless, when imposing Fourier's law in such 
situations, one generally obtains a size-dependent ther- 
mal conductivity k = k(L), which in 1-dimensional sys- 
tems usually in the form of a power-law k oc L a . In addi- 
tion, they observed the absence of a temperature gradient 
inside the system. The same conclusions hold in higher- 
dimensional lattices of classical harmonic oscillators |3J, 
but size-dependence has only been addressed recently, 
while also taking into account disorder and anharmonic- 
ity effects, e.g. in 0HB]. 
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In the present work, however, we focus on the quan- 
tum mechanical equivalent of heat transport though sys- 
tems of harmonic oscillators. While chains of harmonic 
oscillators have been investigated with an emphasis on 
their entanglement properties H], or with the aim 
to scrutinize the domain of validity of master equations 
in Lindblad form to model thermal baths [5], specific 
questions regarding the heat transport properties have 
only been addressed to a limited extend. Within the 
framework of modeling thermal baths by means of quan- 
tum Langevin equations, ballistic transport has been ob- 
served for chains of quantum harmonic oscillators 1 10] . 
More recent analyses include the study of disordered 
chains, but without concrete conclusions regarding the 
exact scaling of the thermal conductivity with the sys- 
tem size Furthermore, explicit analytical results re- 
garding the quantum mechanical steady state remain to 
be elucidated [T5]. We approach this open problem by 
providing analytical solutions to the heat current and an- 
alytical forms of second moments of the non-equilibrium 
steady state in arbitrary dimensions for systems of cou- 
pled harmonic oscillators within the rotating wave ap- 
proximation. 



The paper is structured as follows. In the next section, 
we introduce the employed quantum mechanical model of 
a system of coupled harmonic oscillators, and the master 
equations describing the heat baths. We then proceed in 
section |III| to the analytical solution of the heat current 
in the steady state for a one-dimensional chain. After 
commenting on the thermal nature of the state of indi- 
vidual oscillators in the chain in section |IV| we provide 
the analytical solution for the case that additional local 
dephasing influences the non-equilibrium steady state in 
section [V] before generalizing the result to lattices of ar- 
bitrary dimension in section |Vl| The results are discussed 
in section |VII[ and we summarize in section |VIII| 
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FIG. 1. A linear array of coupled harmonic oscillators with 
thermal environments attached to the endpoints. 

II. MODEL 

The quantum system of interest is a lattice of coupled 
harmonic oscillators. Starting with a chain of N quan- 
tum optical cavities where next-neighbors can coherently 
exchange photons as depicted in Fig. [I] or with an array 
of N pendula coupled by springs, one arrives after the 
rotating-wave approximation at an Hamiltonian of the 
following generic form: 

N N-l 

H = J2 ua ] a 3 + V i^+ 1 ( a K+i + a ]+i a i) > ( 2 ) 
j=i j=i 

where we employ units of h = 1 throughout the paper, w 
denotes the frequency of each of the N identical oscilla- 
tors, Vjj-\-i is the coupling strength between neighboring 
sites, and a,j and at are the bosonic annihilation and 
creation operators, respectively. The interaction term 
thus describes the hopping of a single excitation between 
neighboring sites. Within the performed rotating wave 
approximation, terms proportional to ajOj+i and ajat +1 
have been neglected. The Hamiltonian thus preserves the 
number of excitations and hence commutes with the ex- 
citation number operator. The rotating wave approxima- 
tion is valid in the weak-coupling regime of u > Vjj+i, 
which is often met in quantum optical scenarios. 

The dynamics of the system described by its state p 
is treated within the framework of master equations of 
Lindblad form, that is, the coherent transfer dynamics 
inside the systems on one hand, and the incoherent heat 
exchange with the two baths on the other: 

^ = Cp = -i[H, p] + £ip + C N p. (3) 

Terms C\p and CmP describe the effective, local processes 
of the first and last oscillator of the chain interacting 
with its respective heat bath. Within the Born-Markov 
approximation Cjp can be given in Lindblad form, and 
for a thermal bath it is given by 

Cjp =Tj(rij + 1) (ajpat - ^{a\aj,p}j + (4) 

TjUj {a^jpaj - ^KaJ,pH • 

The first term describes dissipation into the bath, i.e., 
decay of excitations into the reservoir via stimulated 
and spontaneous emission, respectively, and the second 
term describes excitation, i.e., energy absorption from 



the reservoir. The quantity nj — l/[exp(uij/kBTj) — 1] 
denotes the mean excitation number in bath j at the res- 
onance frequency at temperature Tj. This master equa- 
tion of Lindblad form captures the interaction with the 
bosonic heat baths even at high temperatures [5]. 

The expression for the mean heat current through the 
chain can be obtained for the steady state by the follow- 
ing argument. Starting with the time-derivative of the 
energy expectation value, which vanishes in the steady 
state, 

|W -*(*!)-■>, (5) 

one obtains, when inserting the master equation ^ , that 
the coherent part vanishes exactly, while the contribu- 
tions of the heat baths are non-zero but cancel each other. 
The latter amount to the positive net heat current coming 
from the hotter bath and the negative net heat current 
exiting to the colder bath, respectively: 

Tr(H£ 1 p + HC N p) = J 1 + J N =Q. (6) 

The net heat current through the system is thus given 
by J = J i = — Jjv- A compact expression of J for the 
present system employs the specific form of the Hamilto- 
nian ([2| and the terms in the master equation ([11) , and 
when evaluated using commutation relation [cij*,aj] = 1 
yields: 

J = r^i^ni - (ajai)) - ~Vi,2 ((aja 2 ) + (4 a i))- ( 7 ) 

The net heat current through the system is thus char- 
acterized by the mean excitation number of oscillator 1, 
which is coupled to the hotter bath, and by the coher- 
ence between this oscillator and its neighbor, oscillator 2. 
The heat current through the chain is therefore given by 
the difference between the mean energy of oscillator 1 
and the equilibrium state of bath 1, and by the coher- 
ences between the two oscillators at the boundary. An 
analogous expression can be derived for Jjv, in which the 
corresponding terms are expressed in terms of the last 
oscillator, which is connected to the colder bath. 

III. ANALYTICAL SOLUTION 

The coherent dynamics and the incoherent excitation 
exchange with the heat baths in ^ are Gaussian pro- 
cesses that preserve the Gaussian character of Gaussian 
quantum states, i.e. states that can be completely de- 
scribed by their first two moments. Assuming that the 
steady state is the unique solution to dp/dt = Cp = 
we therefore know that the steady state is Gaussian. We 
thus make an ansatz for the non-equilibrium steady state 
using the first and second moments only. Therefore, we 
first define the row vector of all bosonic operators of the 
system 

A = (at, as, . . . ,<Zjv), (8) 
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and the matrix C, which captures the second moments, 
with matrix elements 



[C]ij ^ 



= /J, 



(9) 



which is thus related to A by taking all element-wise 
expectation values of the matrix A^A. 

In order to solve for the moments of the steady state, 
it is first necessary to transform the master equation into 
the equation of motion of the moments. Since the master 
equation is a linear differential equation and calculating 
the moments is also linear in the state, we can treat the 
coherent and incoherent part of the master equation sep- 
arately. The coherent contribution to the equation of 
motion for C can be obtained from the Heisenberg equa- 
tion of motion for A. From the time evolution of the 
annihilation operators, 



in ; 



[aj,H] = Wjdj + Vj,j +l a j+ i + \ , ,. (10) 



in which we formally set a_i = = ajv+i, one obtains a 
system of coupled differential equations, the coefficients 
of which are contained in the matrix W. We thus arrive 
at 



i f = iWA\ 



(11) 



and further, by using the product rule, at the evolution 
equation of A^ A. The coherent part of the differential 
equation for C is thus given by 



dC 

~dt 



i[W,C\. 



(12) 



The effect of the baths can be calculated in a similar way 
by employing the incoherent part of the master equation 
with the exact form of Cj and the bosonic commutation 
relation. We first describe a slightly more general case 
where every oscillator is coupled to a local heat bath via 



d(flj) 
dt 

d(a\aj) 
dt 



d~t~ 



r 

Tr(aiCip) = — ^-(di, 
TT(a\d j [C l p + Cjp}) 



T^aUj) +T j n j 



(13) 



(14) 



(15) 



The incoherent part of the evolution equation of C is thus 
given by 



dC 
~dt 



{L, C} + M, 



(16) 



incoh. 



where the anti-commutator is denoted by {L, C} — LC + 
CL, and the following diagonal matrices were introduced: 



L = -iDia g (r 1 ,r 2) ...,rv), 

M =Diag(rxni,r 2 n 2 , . . .,T N n N ). 



(17) 
(18) 



The complete evolution equation of C for the master 
equation (|3| is thus given by the system of linear dif- 
ferential equations 

dC 

— = i\W,C\ + {L,C} + M, (19) 

where in our scenario of two heat baths at the boundaries 
T 2 = • • • = IV- 1 = 0. Regarding the analytic expression 
of the heat current ^ , the steady state solution of C con- 
tains all the information needed, and all the information 
necessary to describe the Gaussian steady state. 

The ansatz to the solution of C in the steady state, 
dC ss /dt = 0, consists of two contributions: one captur- 
ing the steady state of the average temperature, and one 
that contains the non-perturbative deviations due to the 
non-equilibrium. When there is no temperature bias ap- 
plied, i.e. ri\ — tin = n, the equilibrium steady-state 
solution is given by a diagonal matrix C = nl, i.e. with 
elements [C]u — n on the diagonal and [C]i.j^i = for 
all off-diagonal elements. Each harmonic oscillator of the 
chain thus carries the same average number of excitations 
as both of the heat baths. This case serves as the moti- 
vation for the contribution due to the average excitation 
number. In the non-equilibrium scenario, we thus choose 
the following ansatz 



C ss = nt+ AnD, 



(20) 



where we have defined the average bath excitation num- 
ber n = (n\ +7ijv)/2, and the amount of non-equilibrium 
An = (ni — tin)/2. The matrix D captures the non- 
equilibrium contribution to the steady-state. Insertion 
of this ansatz into (19 1 yields the following equation for 



the non-equilibrium contribution of the steady-state: 

i[D,W}=S + {L,D}, (21) 

with diagonal matrix S = (2nL + M)/An = 
Diag(ri, 0, . . . , 0, Tjv). Given the matrices W with uni- 
form couplings Vjj+i = V, L, and S, the structure of 
(211 implies the following form for the hcrmitian ma- 



trix D, 



D 



x 



V 



x 
e 2 



x 



X 

e N 



(22) 



Thereby, we obtain a set of algebraic equations, which 
for N > 2 reads as follows: 
. V 



i — {x - x*) 
1 1 

»^(ei - e 2 ) 



V 



■ N 



( e j " 
(ejv-i 
. V 



1 - ei, 
x 

~2' 



e j+ i) 
- ejv) 



for 2 < j < N - 2, 



x 
2' 



N 



(X- 



= 1 + e N . 
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The only relevant parameters involved in the solution for 
D are the fractions V/Ti and V/Yn- From the second 
and the second last equation, one can conclude that x is 
purely imaginary, since the e$ are real numbers due to 
the hermiticity of C. Furthermore, with the exception of 
ei and ejv, all the ej are equal. The unique solution to 
the above set of equations reads 



ei 



e 2<j<N-l — 



e N 





4VT 1 r JV 


1 (W 2 + 


rir A r)(r 1 + r A r)' 


4V 2 (T 1 - 


r A r) + r 1 r 2 v + r 2 r JV 


(41/2. 


i-riiv)(ri + iv) 


W 2 (T 1 - 


Tw) + FxT N — V\Y N 


(W 2 - 




W 2 (T 1 - 


Fat) — Tirlr — r 2 rAr 


(W 2 - 


hrxr^^x + rjv) 



The average excitation numbers of the individual oscilla- 
tors on the diagonal of C is thus given by 



{a\ai) 



ei An 



(aja,-) = n + e,, An, 
(ajyaiv) = n + e^An, 



2 < j < N - 1 



that is, the mean excitation number of all oscillators in 
the bulk are equal, while the excitation number of the 
boundary oscillators deviates from that of the bulk to- 
wards that of the heat bath. The coherences between 
neighboring oscillators in the diagonal above and below 
the main diagonal of C are 



(ajoj+i) = xAr 



(23) 



We can now give the analytic expression for the quan- 
tum heat current using the matrix elements of C ss : 



J 



(4V 2 + r 1 r JV )(r 1 + r 



N 



(24) 



The heat current of the chain of quantum harmonic oscil- 
lators is thus independent of the chain length N, and thus 
constitutes a violation of Fourier's law, with the thermal 
conductivity scaling as k ~ N. 

As an alternative expression to Q, we can give the 
heat flux in terms of the purely imaginary coherences, 



J 



2u)Vi(a]jaj 



(25) 



This means, that for the chain of oscillators of the same 
frequency, the heat current is solely given by the next- 
neighbor coherences. We thus expect any additional 
noise source that degrades coherent properties of the 
steady state to degrade the transport properties and de- 
crease the heat current. 



IV. LOCAL THERMALIZATION 

It would be interesting to show that in steady state 
the local density matrices pj at each site are given by a 



thermal state of the following form when written in the 
Fock basis, 



Pj = 



m— 



\m)(m\, (26) 



such that one can locally define a temperature for each 
oscillator in the chain. 

Since the non-equilibrium steady state is of Gaussian 
form, the density operator can be reconstructed from 
the first and second moments. Each individual har- 
monic oscillator of the chain is thus also in a Gaussian 
state. A general expression for the Gaussian density op- 
erator of a single harmonic oscillator [13] is given by 
p = D(a)S(r)pthS(ry D(ay , meaning that any Gaus- 
sian state of this harmonic oscillator can be constructed 
from a thermal state pth by first squeezing it with 



S{r,(j>) = exp(ire~ i2 ^a 2 



and then displac- 



ing it away from the phase space origin with D(a) = 
exp(aa^ — a* a). The mean excitation number of such a 
state is then given by (n) — \a\ +n t h + (2n^ + 1) sinh 2 r, 
where n t h is the thermal contribution to the mean exci- 
tation number. Since due to (13 1 in the steady state 
(cij) = (at) = 0, we expect a to vanish. Furthermore, 
the master equation does not involve any quadratic forms 
in a nor a 1 , and thus we also expect that the squeezing 
operation is not needed to describe the steady state of a 
single oscillator in the chain, i.e. the squeezing parameter 
r should be zero. Under these conditions it is plausible 
that the density operator of a single oscillator at position 
j in the chain is given by a thermal state of Gibbs form 
as was written above. 



TRANSITION BETWEEN COHERENT AND 
INCOHERENT TRANSPORT 



The ballistic transport observed for the uniform non- 
equilibrium chain of quantum harmonic oscillators can 
be turn into diffusive transport obeying Fourier's law by 
adding additional noise terms that degrade next-neighbor 
coherences. Within our master equation approach such 
local noise sources can be straightforwardly implemented 
by introducing additional dephasing terms. Local de- 
phasing environments randomize coherences and thereby 
effectively degrade their magnitude. Thereby, coherently 
delocalized excitations are randomly localized at the sites 
of the harmonic oscillators, which adds a diffusive ele- 
ment to the transport dynamics. The Lindblad super- 
operator for local dephasing processes acting on each of 
the oscillators j is given by 



£*dephP 



N 
3=1 



a/jdjpa^aj 



\ {(a]%) 2 ,p} 



(27) 
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The time evolution of the expectation value of the 
bosonic operators under dephasing is therefore 



dt 

dt 

d(a\ai) 
di 



' 2 



7i 



> i*3, 



0. 



With these equations, the corresponding evolution equa- 
tion of C in the presence of the additional dephasing 
environments gains additional terms: 



dC 
~dt 



dcph 



— {Ldepln C} 

Diag(7i[C]i,i,72[C] 2 ,2, • ■ ■ ,Jn[C]nm) 



where L deph = -|Diag(7i, 72, • • -,7iv)- Let us now con- 
sider the special case in which the local dephasing oper- 
ations acting on the individual oscillators all have equal 
rates, i.e. jj — 7. This allows for a more straightforward 
analytical solution including the dephasing effect. With 
the same ansatz as before, the non-equilibrium part of 
the steady-state solution satisfies the equation 



i[D,W] =S + {L,D}-"/D 

+ 7Diag([ J D] 1>1) [ J D] 2>2) . 



[D}n,n)- (28) 



The algebraic set of equations for the matrix elements 
of D take a similar form as before, but now include the 
dephasing rates: 



V 

i — (x - x*) 
l 1 

« — (ei - e 2 ) 



1 - ei, 



(29) 



i{e 3 



x 
2 



e j+ i) 
- ejv) 



V 



x. 



2 < j < N - 2 



. V 

t N 

v 

i — (x-x*) 

1 N 



X 

~2 ~ I 
1 + e N . 



N 



The solution to the above set of equations are given by 



ei 



t -(4y 2 +r 1 r N )(r 1 +r JV )+2(Jv-i) 7 r 1 r J v ' 

4y 2 (r 1 -r w )+r 1 r 2 v +r 2 r JV +2(iv-i) 7 r 1 r A 
(av 2 +t 1 t n )(t 1 +t n )+2(n-i) 1 t 1 t n 



_ 4\/ 2 (r 1 -r ] y)+r 1 r^-rjr w +2(Ar-2j+i) 7 r 1 r JV 
e-2<j<N-i (4y 2 +r 1 r„)(r 1 +r N )+2(A r -i)7rir N 



e N 



4y 2 (r 1 -r ] y)^r 1 r^-r;r JV -2(jv~i) 7 r 1 r w 
(4y 2 +r 1 r N )(r 1 +r JV )+2(Jv-i) 7 r 1 r JV ■ 



As before, this yields the unique solution to the steady 
state fully contained in C ss . 

In order to derive the analytical expression of the heat 
current in the presence of local dephasing environment, 
it is important to account for a possible heat current 



10" 



-•,=0 
-V=0.01 

•,=0.1 
■ ' 7=0.5 

10 1 



N 




FIG. 2. (Color online.) Left: Heat current vs. chain length 
in a log-log plot for different values of the dephasing rate. 
Right: Mean excitation number Ej = (ajdj) of all oscilla- 
tors for a chain of length iV = 20 in the steady state under 
dephasing. Parameter values: I\2 = 0.1 , V = 0.1, ui — 10. 



due to the dephasing environments. It is a priory not 
clear, that local dephasing does not introduce an addi- 
tional net heat current because the local dephasing pro- 
cesses do not leave energy eigenstates of the chain invari- 
ant. For a chain of harmonic oscillators with uniform 
frequencies, and equal dephasing rates, the net heat cur- 
rent due to the dephasing environment, Tv(HCdephp) = 



■f{a{a 2 ) 



^(a 2 ai) vanishes for the steady state in 



the present scenario because coherences are purely imag- 
inary. Thus the equality J = J\ = — Jjq holds under the 
considered dephasing operations. 

With the non-equilibrium steady state solution for C ss 
the analytical expression for the heat current in the pres- 
ence of dephasing environment is given by 



J = 



4wF 2 r 1 r JV (ni -n N ) 



(w 2 + riivxri + r N ) + 2{n - i) 7 r 1 r 



N 



(30) 



The heat current with local dephasing noise on each os- 
cillator thus acquires a dependence on the size of the 
system. Therefore, the heat current scales as J ~ iV -1 
in the limit of large N, and thereby recovers a size de- 
pendence as in Fourier's law. As before, next-neighbor 
coherences are purely imaginary, but smaller in magni- 
tude. The expression of the heat current in terms of 



coherence (25 1 is also valid in the present case. 



In contrast to the transport scenario without dephas- 
ing, the mean excitations of the individual oscillators in 
the bulk are no longer equal, but decrease linearly from 
the hotter towards the colder heat bath. The effects of 
dephasing on the individual mean excitations of the os- 
cillators is shown in Fig. [2] In the limit of large dephas- 
ing rates, the entire harmonic chain approaches the phe- 
nomenology of Fourier's law, in the sense that a constant 
gradient emerges throughout the chain. 



VI. GENERALIZATION TO d-DIMENSIONAL 
LATTICES 

With the analytical results for heat transport through 
chains (dimension d — 1) and the strategy to solve for 
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FIG. 3. A 2-dimensional harmonic lattice whose boundary 
oscillators of opposite edges (left and right) are coupled to 
local heat baths of different temperature. 



steady state properties available, we generalize our re- 
sults to heat transport scenarios through lattices of d 
dimensions. The Hamiltonian ^ is extended to a cu- 
bic lattice of d dimensions, again with uniform oscilla- 
tor frequencies for all oscillators, and uniform couplings 
between next neighbors along all lattice edges. When 
advancing to higher dimensions, we consider a situation 
where the heat transport is only driven along one direc- 
tion. Therefore, two opposing (d — l)-dimensional hyper- 
surfaces are coupled to heat baths of different tempera- 
ture such that there is a temperature difference along the 
one remaining dimension. The two heat baths are mod- 
eled such that all harmonic oscillators of the respective 
hyper-surface are coupled to a local heat bath of iden- 
tical temperatures. An example, for the d = 2 case is 
depicted in Fig. [3] with the two opposing edges coupled 
to local heat baths. 

We approach the analytic solution by relating it to the 
1-dimensional case using a recursion argument. Clearly, 
due to symmetry we do not expect a net heat current 
perpendicular to the applied temperature gradient. How- 
ever, the coherent coupling in lattice directions trans- 
verse to the temperature gradient might constructively 
or destructively interfere with the heat transport in the 
longitudinal direction. As for the 1-dimensional case, 
we assume the uniqueness of the non-equilibrium steady 
state. This again leads us to the solution, with all es- 
sential information being captured in a matrix C with 
matrix elements (a\a,j), where i and j are index vec- 
tors that account for labeling of all harmonic oscillators 
in d dimensions. We choose the order of indices such 
that the last index counts oscillators along the direction 
along which the temperature bias is applied. Thereby we 
can relate the heat transport through a d-dimensional 
lattice to the heat transport through d stacked copies 
of (d — l)-dimensional lattices that are coherently cou- 
pled. The matrix C thus is of block structure. For a d- 
dimensional hypercube of N sites along every direction, 
C = Cd can be divided into N x N blocks, where each 
block of the diagonal captures a matrix G d _ x of the j'-th 
(d — l)-dimensional lattice, and the off-diagonal blocks 
capture the coherences between these different lattices 
of lower dimension. Since the stacked lower-dimensional 
lattices are only coupled to their direct neighbors, only 



coherences between two neighboring lattices will exists, 
hence only the blocks on the first diagonals above and be- 
low the main diagonal being non-zero. This procedure is 
repeated until the two-dimensional lattice in decomposed 
into an array of coupled chains. 

Similarly, all the sets of differential equations can be 
related to those of the next-lower dimension with ad- 
ditional couplings along the transverse direction. We 
choose the same ansatz for the non-equilibrium steady 
state, Cd — n\ + AnDd, which when inserted into the 
master equation yields 

i[D d ,W d ]=S d + {L d ,D d }. (31) 

Let us now employ the described block structure 



(D d % Q 



D ri 



\ 



V 



Q* d£?> q 



(Wd-i g 
G Wd-i G 



Wd 



\ 



V 



G W d -i G 
G W d -J 



where G — VTL is a diagonal matrix capturing the cou- 
pling of harmonic oscillators between two neighboring 
(d— l)-dimensional sub-lattices. The matrix Q accounts 
for the coherences between these oscillators. The matri- 
ces S d and Ld again capture the influence of the heat 
bath, and they are simply given by 

N N 

L d = ($L d - 1 and S d = 0S d _i, (32) 

3=1 3=1 

with L\ and Si being those for the linear chain given in 

sect, nn 



Inserting these matrices into (311 gives rise to the set 



of differential equation for the lower-dimensional sub- 
lattices: 

i[D® v Wd-!} = S d -! + {Ld-LD^} Vj, (33) 
{L d -i,Q d } = 0. (34) 

These equations imply that Qd vanishes and hence there 
are no coherences and no heat transport between sub- 
lattices. Furthermore, since the equations are identical 
for all D d _ 1} the steady state is given by the identical 
steady states of the lower-dimensional sub-lattice. The 
non-equilibrium steady state of the d-dimensional lat- 
tice is thus given by the steady states of the chains con- 
necting heat baths at different temperature to which we 
have decomposed the lattice. The heat current of the 
d-dimcnsional lattice is thus given by the sum of all the 
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heat currents of the identical chains composing the lat- 
tice. Hence the total heat current in tf-dimension is the 
volume of the hyper-surfaces that are in contact with the 
heat baths, i.e. the number of oscillators coupled to heat 
baths of the same temperature, times the heat current of 
the chain as derived in Sec. |III| For the special case of a 
d-dimcnsional hypercube with an edge length of N sites 
in each direction 



J d = Yol d ^J = N d ' 1 J. 



(35) 



This implies that the heat current per chain is invariant 
under a length change in the direction parallel to the heat 
current. The heat current of this model is ballistic also 
in higher dimensions. 



VII. DISCUSSION 

The analytic results for heat transport in harmonic lat- 
tices of arbitrary dimension show a number of differences 
to existing work. The framework of master equations 
in Lindblad form readily provides access to not only to 
steady state properties, but to the steady state itself, 
and it is straightforwardly applicable to many quantum 
optical setups and scenarios, such as the array of cou- 
pled cavities depicted in Fig. [T] Our results differ from 
that of a harmonic chain where the baths are modeled by 
quantum Langevin equations |10L 111] . On one hand the 
method of treating the heat baths differs, on the other 
hand the ballistic transport observed in our model is in- 
dependent of the system size in the absence of dephasing. 
Furthermore, while the mean excitation profile of the os- 
cillators coincides in some parameter regimes, we obtain 
different profiles in other cases. We attribute these dif- 
ferences to the performed rotating wave approximation 
in the Hamiltonian of our model, but a detailed investi- 
gation of the effects of the rotating wave approximation 
on the heat transport is beyond the scope of the present 
work. 

With respect to the open question regarding explicit re- 
sults for the steady state density matrix as stated in [12] , 
our work supplies analytical answers for the special case 
of the considered Hamiltonian and bath models, even in 
arbitrary dimension. The simple relation of the heat cur- 
rent in lattices of higher dimension may serve as a guiding 
result for numerical studies in higher dimensions. 

Approaching the diffusive transport regime of Fourier's 
law by adding dephasing environments is consistent with 
other studies [TH - fT?] . which investigate the conditions 
under which the heat transport in quantum systems ap- 
proaches a diffusive regime obeying Fourier's law. The 
mechanism applied in [ il4] in order to obtain heat dif- 
fusion in quantum chains, i.e. introducing a band of ex- 
cited states rather than sharp energy levels per site, can 
qualitatively be recovered by our dephasing process, i.e., 
the effectively introduced level fluctuations and hence the 
line broadening as also observed for chains of two-level 
atoms |16| . The dynamic disorder introduced by the local 



dephasing environments allows for obtaining analytic re- 
sults and a more straightforward clarification of the size- 
dependence than approaching the diffusive regime with 
disorder in the Hamiltonian as pursued in 

In transport scenarios of molecular biology, e.g. in 
light-harvesting systems of photosynthesis, it is being in- 
vestigated to what extent entanglement or (in general 
non-equivalently) measures of coherence may serve as an 
indicator of transport properties |18If22] . While we have 
observed the emergence of entanglement in strongly non- 
equilibrium systems of coupled two- level atoms |16j . it 
does not hold for systems of harmonic oscillators. 

This result is not surprising given that operators of the 
form a 2 and (a^) 2 have been neglected in the Hamilto- 
nian and the master equation, which are required, how- 
ever, to achieve a squeezing operation that is necessary 
for creating an entangled Gaussian state. The conclu- 
sion to be drawn is that although the system can exhibit 
large coherences and thereby a large heat current, i.e. ef- 
ficient transport, entanglement is completely absent and 
also cannot be generated by a large temperature differ- 
ence as possible for systems of two- level atoms [T3] . As 
spelled out by ( 25 ) , the magnitude of the heat current 



depends on the amount of quantum coherence between 
the neighboring oscillators. Therefore, entanglement cer- 
tainly does not qualify as a signature of quantum trans- 
port efficiency in the present model. 

The validity of Fourier's law in higher dimensions 
strongly depends on the modeling of the system. Most 
of the works in this direction highly rely on large-scale 
simulations [SJ |SJ [23] . Peierls showed that the phonon- 
phonon scattering leads to diffusive heat conduction in 
higher dimensional systems [25]. In [25] it is shown, that 
all models characterized by short range interactions and 
momentum conservation should exhibit the same kind of 
anomalous behavior in the heat conduction, i.e. a diver- 
gence of the thermal conductivity at infinite system-size, 
for d < 3. We can thus extend the result of [23] for the 
specific case of harmonic oscillator lattices to arbitrary 
dimensions d with an anomalous heat conduction with 
k oc N. 



VIII. SUMMARY 

We provide an analytical solution to the heat cur- 
rent and the non-equilibrium steady state of a chain of 
quantum harmonic oscillators, whose boundary oscilla- 
tors are coupled to two heat baths of different temper- 
ature, respectively. The heat current is independent of 
the chain length, and thus is of ballistic nature exhibit- 
ing an anomalous heat current with a thermal conduc- 
tivity proportional to the chain length. The diffusive 
heat transport regime, i.e. normal heat conduction with 
a constant thermal conductivity, is recovered when addi- 
tional dephasing environments locally affect each of the 
oscillators. We observe the absence of entanglement in 
the chain for all parameter regimes. Finally, we pro- 
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vide the analytical expression for the heat current of a 
rf-dimcnsional lattice of harmonic oscillators, which turns 
out to be the sum of the heat currents of all the chains, 
into which the lattice may be decomposed, each one con- 
necting the two heat baths. 
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